Graph Analysis with R. Part B

Description

In this project, we were provided with a raster file (.tif) containing elevation data for the Mediterranean, along with three shape files containing information about capitals, states, and places in Greece. The project involves various geocomputations using these files, and new data will be added to create additional interesting maps by the end.

Add libraries 

library(sf)
library(terra)
library(dplyr)
library(spData)
library(raster)
library (ggplot2)
library(ggrepel)
library(dplyr)
library(stringr)

Read and plot raster and vector files 

#raster data
file_path <- "myrast.tif"
raster_data <- raster(file_path)
print(raster_data)
## class      : RasterLayer 
## dimensions : 2400, 3600, 8640000  (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : -0.0001388889, 29.99986, 29.99986, 49.99986  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : myrast.tif 
## names      : myrast
plot(raster_data, main = "Mediterranean map")

#shape data
poleis = st_read("poleis/poleis.shp")
## Reading layer `poleis' from data source `/users/aid24005/poleis/poleis.shp' using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 149581.4 ymin: 3895170 xmax: 878802.9 ymax: 4555185
## Projected CRS: GGRS87 / Greek Grid
places = st_read("places/places.shp")
## Reading layer `places' from data source `/users/aid24005/places/places.shp' using driver `ESRI Shapefile'
## Simple feature collection with 4326 features and 4 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 19.39397 ymin: 34.84761 xmax: 27.2865 ymax: 41.72253
## Geodetic CRS:  WGS 84
nomoi = st_read("GRC_ADM2/GRC_ADM2.shp")
## Reading layer `GRC_ADM2' from data source `/users/aid24005/GRC_ADM2/GRC_ADM2.shp' using driver `ESRI Shapefile'
## Simple feature collection with 77 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 19.3726 ymin: 34.8009 xmax: 29.60684 ymax: 41.74889
## Geodetic CRS:  WGS 84

Raster Cropping Greece 

###1
places <- st_transform(places, crs(raster_data)) #From Geocomputation with R ,6.2 Raster cropping

greece_geometry <- st_union(places$geometry) #Compute Greece Geometry from union its states

greece_raster <- crop(raster_data, st_bbox(greece_geometry)) #Crop Greece from initial raster

plot(greece_raster, main= "Greece map")

Elevation of states capitals

A map of the regions of Greece, where the size of each point corresponding to each capital is proportional to the elevation of that capital. The graph on the left represents each state’s capital, and on the right, the size of the circle for each capital depends on its elevation

###2

#compute capitals altitude
poleis <- st_transform(poleis, crs(greece_raster))
poleis$altitude<- extract(greece_raster , poleis) #capitals altitude

#create map with states and state's capitals 
par(mfrow = c(1, 2)) #2 plots
cex = sqrt(poleis$altitude)/10
plot(nomoi["NAME"] , reset = FALSE, main = "Capitals of states")
plot(st_geometry(poleis), add = TRUE) 

plot(nomoi["NAME"] , reset = FALSE, ,main="Elevation of capitals") 
plot(st_geometry(poleis), add = TRUE , cex=cex) #The biger the circle is, the maximum altitude it has

par(mfrow = c(1, 1))

Mean elevation for states 

###3
altitude_nomoi<-extract(raster_data, nomoi) #Computes elevations for all the state 
nomoi$mean_altitude <- sapply(altitude_nomoi, mean, na.rm = TRUE) #mean
nomoi$sd_altitude <- sapply(altitude_nomoi, sd, na.rm = TRUE) #standart deviation
plot(nomoi[9], main = "Nomoi - mean_elevation")

Standard Deviation for the elevation of states 

plot(nomoi[10], main= "Nomoi - standard_deviation_elevation")

The absolute value of the difference between the mean elevation of the state and the elevation of its capital

In this step, it is necessary to create a new data frame combining features from both ‘nomoi’ and ‘poleis’. Then, compute the difference between capitals and states and represent them in an improved plot.

###4
# Fix invalid geometries
nomoi <- st_make_valid(nomoi)
poleis <- st_make_valid(poleis)

#intersection of poleis, nomoi
intersection_nomoi_capitals <- st_intersection(poleis, nomoi) 

#join it with poleis to add poleis$altitude
nomoi_capitals_pairs <- st_join(intersection_nomoi_capitals, poleis) 

# Calculate the absolute difference
merged <- nomoi_capitals_pairs %>%
  mutate(difference = abs(mean_altitude- altitude.x))

#Convert as sf object
sf_merged <- st_as_sf(merged)

#join with nomoi, so it is easier to plot colored nomoi geometry
plot_nomoi <- st_join(nomoi, merged)

# Plot choropleth map
ggplot(plot_nomoi) +
  geom_sf(aes(fill = difference)) +
  scale_fill_gradient(low = "blue", high = "red", name = "Absolute Altitude Difference") +
  ggtitle("Choropleth Map of Absolute Altitude Difference") +
  theme_minimal()

The top 10 states in terms of Average Elevation 

###5a
#filter for top 10 states by mean_altitude
top_10_nomoi_mean_altitude <- nomoi %>% 
  arrange(desc(mean_altitude)) %>%  
  head(10)  

print(paste("Top 10 states by its mean altitude are:" ))
## [1] "Top 10 states by its mean altitude are:"
print(top_10_nomoi_mean_altitude$NAME)
##  [1] "Kastoria Regional Unit"     "Regional Unit of Evrytania"
##  [3] "Florina Regional Unit"      "Ioannina Regional Unit"    
##  [5] "Grevena Regional Unit"      "Regional Unit of Phocis"   
##  [7] "Regional Unit of Kozani"    "Arcadia Regional Unit"     
##  [9] "Regional Unit of Trikala"   "Drama Regional Unit"

The top 10 states in terms of Standard Deviation of Elevation 

###5b
#filter top 10 states by standard deviation
top_10_nomoi_sd_altitude <- nomoi %>% 
  arrange(desc(sd_altitude)) %>%  
  head(10)  

print(paste("Top 10 states by its mean altitude are:" ))
## [1] "Top 10 states by its mean altitude are:"
print(top_10_nomoi_sd_altitude$NAME)
##  [1] "Imathia Regional Unit"    "Pieria Regional Unit"    
##  [3] "Regional Unit of Trikala" "Pella Regional Unit"     
##  [5] "Regional Unit of Phocis"  "Chania Regional Unit"    
##  [7] "Arta Regional Unit"       "Achaea Regional Unit"    
##  [9] "Karditsa Regional Unit"   "Ioannina Regional Unit"

Places that are situated above 1500 meters in altitude 

###6
places$altitude<- extract(greece_raster , places) #places altitude

#filter for places with altitude>1500
places_over_1500 <-places %>%
  filter(altitude>1500)

#initial plot of all states
base_plot <- ggplot() +
  geom_sf(data = nomoi[0]) 

#plot by population difference
second_plot<-base_plot+
  #add places over 1500, colored by its population
  geom_sf(data = places_over_1500, aes(color =population > 0))+
  #Process the text labels to avoid overlap and ensure a visually pleasing appearance
   geom_text_repel(
    data = places_over_1500,
    aes(x = st_coordinates(geometry)[, 1], y = st_coordinates(geometry)[, 2], label = name),
    box.padding = 0.2, 
    point.padding = 0.1,  
    segment.color = NA,
    segment.size = 0,
    direction = "both",  
    hjust = 0, vjust = 0
  )


plot(second_plot)

Reclassify Greece map 

###7
rcl <- matrix(c(1, 500, 1, 500, 1000, 2, 1000, 1500, 3, 1500, 2000, 4, 2000, 2500, 5, 2500, 3000, 6), ncol = 2 , byrow = TRUE) #create reclassifies levels
recl <- reclassify(greece_raster, rcl = rcl) #reclassify map
plot(recl, main= "Reclassified Greece map") #plot

Line between Grebena-Ioannina and elevation along the line 

###8
grevena<-poleis %>%
  filter(NAME == "Grebena")
ioannina<-poleis %>%
  filter(NAME =="Ioannina")

# Extract coordinates for Grevena
coordinates_grevena <- st_coordinates(grevena$geometry)
x_value_grevena <- coordinates_grevena[, "X"]
y_value_grevena <- coordinates_grevena[, "Y"]

# Extract coordinates for Ioannina
coordinates_ioannina <- st_coordinates(ioannina$geometry)
x_value_ioannina <- coordinates_ioannina[, "X"]
y_value_ioannina <- coordinates_ioannina[, "Y"]

# Δημιουργία γραμμής που ενώνει τις δύο πόλεις
greece_transect1 <- st_linestring(rbind(c(x_value_grevena, y_value_grevena) , c(x_value_ioannina, y_value_ioannina))) %>%
  st_sfc(crs = st_crs(greece_raster)) |>
  st_sf(geometry = _)


#Geocomputation with R: 5.4.2 Raster extraction

greece_transect1$id = 1:nrow(greece_transect1)
greece_transect1 = st_segmentize(greece_transect1, dfMaxLength = 50)
greece_transect1= st_cast(greece_transect1, "POINT")

greece_transect1 = greece_transect1 |> 
  group_by(id) |> 
  mutate(dist = st_distance(geometry)[, 1]) 

greece_transect1_elev = terra::extract(greece_raster, greece_transect1)
greece_transect1 = cbind(greece_transect1, greece_transect1_elev)


greece_transect1_elevation_data <- data.frame(elevation = greece_transect1_elev)

par(mfrow = c(1, 2), mar = c(1, 6, 5, 2)) #2 plots, avoid overlaping between them

#1
plot(greece_raster , main= "Line extraction: Ioannina - Grebena")
plot(greece_transect1, add = TRUE, col = "red")

#2
plot(greece_transect1_elevation_data$elevation, type = "l", col = "blue", xlab = "Distance", ylab = "Elevation", main= "Elevation along the line")

par(mfrow = c(1, 1))
###8
katerini<-poleis %>%
  filter(NAME == "Katerini")
lamia<-poleis %>%
  filter(NAME =="Lamia")

# Extract coordinates for Katerini
coordinates_katerini <- st_coordinates(katerini$geometry)
x_value_katerini <- coordinates_katerini[, "X"]
y_value_katerini <- coordinates_katerini[, "Y"]

# Extract coordinates for Lamia
coordinates_lamia <- st_coordinates(lamia$geometry)
x_value_lamia <- coordinates_lamia[, "X"]
y_value_lamia <- coordinates_lamia[, "Y"]

# Δημιουργία γραμμής που ενώνει τις δύο πόλεις
greece_transect2 <- st_linestring(rbind(c(x_value_katerini, y_value_katerini) , c(x_value_lamia, y_value_lamia))) %>%
  st_sfc(crs = st_crs(greece_raster)) |>
  st_sf(geometry = _)


#Geocomputation with R: 5.4.2 Raster extraction

greece_transect2$id = 1:nrow(greece_transect2)
greece_transect2 = st_segmentize(greece_transect2, dfMaxLength = 50)
greece_transect2= st_cast(greece_transect2, "POINT")

greece_transect2 = greece_transect2 |> 
  group_by(id) |> 
  mutate(dist = st_distance(geometry)[, 1]) 

greece_transect2_elev = terra::extract(greece_raster, greece_transect2)
greece_transect2 = cbind(greece_transect2, greece_transect2_elev)


greece_transect2_elevation_data <- data.frame(elevation = greece_transect2_elev)

par(mfrow = c(1, 2), mar = c(1, 6, 5, 2)) #2 plots, avoid overlaping between them

#1
plot(greece_raster , main= "Line extraction: Katerini - Lamia")
plot(greece_transect2, add = TRUE, col = "red")

#2
plot(greece_transect2_elevation_data$elevation, type = "l", col = "blue", xlab = "Distance", ylab = "Elevation", main= "Elevation along the line")

par(mfrow = c(1, 1))

Playgrounds in Athens and their distances from the city center

For the creation of the plot, a shape file from the link was used, containing the playgrounds in the Athens region. By utilizing the st_centroid function for Central Athens, the visual representation shows the distance of these playgrounds from the city center.

###9a
nomos_athens <- nomoi %>%
  filter(str_detect(NAME, "Central Athens"))

paidikes_xares <- st_read("paidikes_xares_da/paidikes_xares_da.shp")
## Reading layer `paidikes_xares_da' from data source 
##   `/users/aid24005/paidikes_xares_da/paidikes_xares_da.shp' using driver `ESRI Shapefile'
## Simple feature collection with 28 features and 6 fields
## Geometry type: MULTIPOINT
## Dimension:     XY
## Bounding box:  xmin: 474045.8 ymin: 4200701 xmax: 479752.4 ymax: 4208680
## Projected CRS: GGRS87 / Greek Grid
# http://geodata.gov.gr/en/dataset/tt1ototto1nueves-tta1d1kes-xapes-anuou-a0nvaiwv


athens_centroid <- st_centroid(nomos_athens)


# Σχεδίαση του χάρτη
ggplot() +
  geom_sf(data = nomos_athens, fill = "lightblue", color = "black") +
  geom_sf(data = paidikes_xares, fill = "red", color = "red") +
  geom_sf(data = athens_centroid, color = "blue", size = 3) +
  ggtitle("Playgrounds Distance to Athens Centroids")

Railroads in Greece

Here, a shapefile containing railroads of Greece was employed. By following the provided link and selecting Greece and Railroads. The first graph displays the railroads in Greece. In the subsequent graph, the states of Greece are represented with colors corresponding to the count of different railroads they have.

###9b

railroads <- st_read("GRC_rails.shp") #https://www.diva-gis.org/datadown
## Reading layer `GRC_rails' from data source `/users/aid24005/GRC_rails.shp' using driver `ESRI Shapefile'
## Simple feature collection with 133 features and 7 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 21.11253 ymin: 37.04311 xmax: 26.63872 ymax: 41.73741
## Geodetic CRS:  WGS 84
railroads_nomoi_intersection <- st_intersection(railroads, nomoi) #intersection of which nomoi have railroads

# Convert sf objects to data frames 
railroads_df <- st_as_sf(as.data.frame(railroads))
nomoi_df <- st_as_sf(as.data.frame(nomoi))
railroads_nomoi_intersection_df <- st_as_sf(as.data.frame(railroads_nomoi_intersection))

# Count times of nomoi that appears in Dataset interse
nomoi_counts <- railroads_nomoi_intersection_df %>%
  group_by(NAME) %>%
  summarise(count = n())

nomoi_counts_df <- st_drop_geometry(nomoi_counts)

# Use st_join to add the counts to the nomoi_df
nomoi_df <- left_join(nomoi_df, nomoi_counts_df, by = "NAME")


ggplot()+
  geom_sf(data = nomoi_df) +
  geom_sf(data =railroads ) +
  ggtitle("Railroads in Greece")

ggplot() +
  geom_sf(data = nomoi_df, aes(fill = count), color = "black") +
  scale_fill_gradient(low = "green", high = "red") +
  ggtitle("Number of Railroads Passing Through Each State")

Roads of Greece

In this section, the roads of Greece are represented using the roads link. By selecting Greece and Roads in the provided link, the roads consist of linestrings and multilinestrings. Using st_length, the longest road in Greece is determined and plotted.

###9c

roads = st_read("GRC_roads.shp") #https://www.diva-gis.org/datadown
## Reading layer `GRC_roads' from data source `/users/aid24005/GRC_roads.shp' using driver `ESRI Shapefile'
## Simple feature collection with 2222 features and 5 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 19.64183 ymin: 34.82458 xmax: 28.22389 ymax: 41.73907
## Geodetic CRS:  WGS 84
#find roads length
roads$length<- st_length(roads$geometry)
#max length
max_length_index <- which.max(roads$length)

# Extract the feature with the maximum length
max_road <- roads[max_length_index, ]

#Max road in which take place only in one nomo

ggplot()+
  geom_sf(data = nomoi_df) +
  geom_sf(data =max_road) +
  ggtitle("Longest road")

As observed, consider the maximum road length, is in Athos, where we currently lack additional data from the nomoi.shp file. Therefore, we will utilize the second-longest road.

#second maximum road
second_max_length_index <- which.max(roads$length[-max_length_index])

# Extract the feature with the second maximum length
second_max_road <- roads[second_max_length_index, ]

ggplot()+
  geom_sf(data = nomoi_df) +
  geom_sf(data =second_max_road) +
  ggtitle("Second longest road ")

print(second_max_road)
## Simple feature collection with 1 feature and 6 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 22.01608 ymin: 40.1975 xmax: 22.32342 ymax: 40.40264
## Geodetic CRS:  WGS 84
##     MED_DESCRI RTT_DESCRI F_CODE_DES ISO ISOCOUNTRY
## 371       <NA>       <NA>      Trail GRC     GREECE
##                           geometry       length
## 371 MULTILINESTRING ((22.32342 ... 66736.34 [m]

Using st_intersection, the states to which this road belongs are color-coded in the plot

intersection_road_nomoi <- st_intersection(second_max_road, nomoi)

selected_nomoi <- intersection_road_nomoi$NAME
selected_nomoi_data <- nomoi %>%
  filter(NAME %in% selected_nomoi)

ggplot()+
  geom_sf(data = nomoi_df) +
  geom_sf(data = selected_nomoi_data, map=aes(fill = "blue")) +
  geom_sf(data =second_max_road, map = aes(fill= "black")) +
  ggtitle("The state in which the second longest road is located.") +
  theme_minimal() +
  theme(
    legend.position = "none")

Usefull links Geocomputation with R